/*!
  \file
  Implementation of the FFT for HHVolumeSimple (uses FFTW)
*/

#include <fftw3.h>

#include "fftimagesimple.hpp"


// =========================================

void FFTimageSimple::Forward( SciGPU::HHV::HHVolume< std::complex<float>,2>& dst,
                              const SciGPU::HHV::HHVolume<float,2>& src ) const {
  SciGPU::HHV::HHVolumeSimple< float,2>::convertor srcConvert( src );
  SciGPU::HHV::HHVolumeSimple< std::complex<float>,2>::convertor dstConvert( dst, false );
  
  this->ForwardTransform( dstConvert.get(), srcConvert.getConst() );
}




void FFTimageSimple::ForwardTransform( SciGPU::HHV::HHVolumeSimple< std::complex<float>,2>& dst,
                                       const SciGPU::HHV::HHVolumeSimple<float,2>& src ) const {
  /*!
    This does a very lazy transform - convert the input to complex and do a full transform
  */
  SciGPU::HHV::HHVolumeSimple< std::complex<float>,2> srcComplex;
  srcComplex.Convert( src );

  
  const SciGPU::HHV::DimensionObject<2> dims = src.dimensions();


  dst.Allocate( dims );

  fftwf_plan thePlan;
  thePlan = fftwf_plan_dft_2d( dims[0], dims[1],
                               reinterpret_cast<fftwf_complex*>(srcComplex.start()),
                               reinterpret_cast<fftwf_complex*>(dst.start()),
                               FFTW_FORWARD,
                               FFTW_ESTIMATE|FFTW_PRESERVE_INPUT|FFTW_UNALIGNED );
  if( NULL == thePlan ) {
    std::cerr << __FUNCTION__
              << ": Failed to create plan"
              << std::endl;
    abort();
  }

  // Run it
  fftwf_execute( thePlan );

  // Release the plan
  fftwf_destroy_plan( thePlan );
  
}



// ==========================================


void FFTimageSimple::Inverse( SciGPU::HHV::HHVolume<float,2>& dst,
                              const SciGPU::HHV::HHVolume< std::complex<float>,2>& src ) const {
  SciGPU::HHV::HHVolumeSimple< float,2>::convertor dstConvert( dst, false );
  SciGPU::HHV::HHVolumeSimple< std::complex<float>,2>::convertor srcConvert( src );
  
  this->InverseTransform( dstConvert.get(), srcConvert.getConst() );
}


void FFTimageSimple::InverseTransform( SciGPU::HHV::HHVolumeSimple<float,2>& dst,
                                       const SciGPU::HHV::HHVolumeSimple< std::complex<float>,2>& src ) const {
  const SciGPU::HHV::DimensionObject<2> dims = src.dimensions();
  SciGPU::HHV::HHVolumeSimple< std::complex<float>,2> dstComplex;

  dstComplex.Allocate( dims );
  dst.Allocate( dims );


  fftwf_plan thePlan;

  thePlan = fftwf_plan_dft_2d( dims[0], dims[1],
                               reinterpret_cast<fftwf_complex*>(const_cast<std::complex<float>*>(src.start())),
                               reinterpret_cast<fftwf_complex*>(dstComplex.start()),
                               FFTW_BACKWARD,
                               FFTW_ESTIMATE|FFTW_UNALIGNED|FFTW_PRESERVE_INPUT );
  if( NULL == thePlan ) {
    std::cerr << __FUNCTION__
              << ": Failed to create plan"
              << std::endl;
    abort();
  }


  // Run it
  fftwf_execute( thePlan );
  
  // Release it
  fftwf_destroy_plan( thePlan );

  // Rescale and convert data
  const std::complex<float> * const dstCpxArr = dstComplex.start();
  float * const dstArr = dst.start();
  for( unsigned int i=0; i<dstComplex.hyxels(); i++ ) {
    dstArr[i] = dstCpxArr[i].real() / ( dims[0] * dims[1] );
  }
}
  
